##################################################################
##
## Conflict in Lake Chad region
##
## B. Blankespoor 2020-08-16 (first version 2020-04-01)
##
##################################################################


options(rasterTmpDir='C:/temp/')
roads_dir <- 'V:/01_SSA/018_transportation/REMI_JEDWAB_restricted/GISroads_Adam_Remi'
## Slope_globe<-raster("V:/00_GLB/006_elevation/slope/meanslpprc/meanslpprcg")
# #Slope_globe<-raster("E:/wb/data/GLB_00/006_elevation/meanslope/meanslope.tif")
afripolis_shp = 'V:/01_SSA/016_society/Africapolis_2015/africapolis'
  
##
## IMPORTANT PATH to Rtools to run Cpp
##
Sys.setenv(PATH = paste("C://WBG//Utility//Rtools//bin", Sys.getenv("PATH"), sep=";"))
Sys.setenv(BINPREF = "C://WBG//Utility//Rtools//mingw_$(WIN)//bin//")

require(Rcpp)
output_version <- Sys.Date()
  
##
## UDF
##


#Tracks, improved, paved, highways and speed (r) but do Euclidean distance only (d) => two measures will be different if two measures are not strongly correlated. 
#Highway = 80. paved = 60. improved = 40. dirt = 12. no road = 6
#Assumption = free crossing (open border) (o) (later on, we will have p = porous, p1 p2.)
#Sigma = 1, 2, 3, 4 (1 2 3 4)
#For example, baseline = lcmacro4 for each cell-year
#RXXXX: Surface of the road in the year XXXX (e.g., 1986 or 2003).
#There are four surfaces:
#8 = highway (highways are paved roads but with at least 2 x 3 lanes) 
#1 = paved road
#3 = improved road
#5/0 = earthen road


## speeds in kph
## speeds in kph
roadcat2speed <- function(road_cat) {
  
  if (road_cat==8) {
    return(80)
  }
  if (road_cat==1) {
    return(60)
  }
  if (road_cat==3) {
    return(40)
  }
  if (road_cat==5) {
    return(12)
  } 
  if (road_cat==0) {
    return(6)
  } else {
    return(6)
  }
}

calcPoints <- function(acled_loc){
  # add small value when no ntl #
  acled_loc[acled_loc==0] = 1
  rast.df <- as.data.frame(rasterToPoints(acled_loc))  
}


calcConflictAccessParallel <<- function(cost.t,fromCoordsWeight,sizegrid,trade_param,.packages=c("gdistance","raster","igraph", "Matrix")){
  nccell <- Which(sizegrid==1, cells = TRUE)
  sizeconfgrid = rasterFromXYZ(fromCoordsWeight, res=res(sizegrid), crs=crs(sizegrid))
  # match extents to country, conf grid only from points
  sizeconfgrid = extend(sizeconfgrid,sizegrid)
  
  sizegridid = sizegrid
  values(sizegridid) = 1:ncell(sizegrid)
  
  # multiply by 1 and NA to reduce grid #
  sizegridid = sizegridid * sizegrid

  ## SUBSET LOCATIONS ONLY WITH DATA ON GRID  ##
  cum_sum <- calc(sizegridid, fun = function(x) {
               fromCoordsMat = xyFromCell(sizegridid, x, values=FALSE)
               # travel time cost
               if (dim(fromCoordsMat)[1]==1) {
                 # sizeconfgrid - The function uses Dijkstra's algorithm (as implemented in the igraph package).
                 acc_cost.hour <- gdistance::accCost(cost.t, fromCoordsMat)
                 acc_cost.hour[acc_cost.hour==0] <- NA
               # Exclude origin, which is equal to 0
               # Conf [event,fatal,combined] * [night time lights] * (travel time in hours)^(trade_param)
               acc_cost.ma.hour <- (sizeconfgrid * acc_cost.hour)^-trade_param
               cma <- cellStats(acc_cost.ma.hour, stat='sum', na.rm=TRUE)
               return(cma)
               }
             })
  return(cum_sum)
}

adm0_lake_chad_shp = paste0(wkdir,'/local/gis_data/003_boundaries/lake_chad/lake_chad_adm0')
adm0_reg_shp <- st_read(dirname(adm0_lake_chad_shp), basename(adm0_lake_chad_shp))
adm0_lake_chad_sf <- read_sf(dsn=dirname(adm0_lake_chad_shp), layer=basename(adm0_lake_chad_shp))
fishnet_ext <- extent(adm0_reg_shp)

##
## BH between rivers
##
bhbr_shp <- paste0(wkdir,'/local/gis_data/003_boundaries/bhbr_dst')
bhbr <- read_sf(dsn=dirname(bhbr_shp),layer=basename(bhbr_shp))

##
## READ PROCESSED FISHNET HERE
##
## fishnet_only
## select by location intersection of countries
## add dummy
## feature to point
## point to raster (tif)
print('fishnet only - no adm boundaries')
fishnet_lake_chad_ntl_x_d01_shp <- paste0(wkdir,'/local/gis_data/003_boundaries/fishnet/fishnet_lake_chad_ntl_x_d01')
fishnet_only <- read_sf(dirname(fishnet_lake_chad_ntl_x_d01_shp),basename(fishnet_lake_chad_ntl_x_d01_shp))
fishnet_tif <- paste0(wkdir,'/local/gis_data/003_boundaries/fishnet/fishnet_lake_chad_ntl_x_d01.tif')
fishnet_grid <- raster(fishnet_tif)

ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01')
names(ntl.zones)
summary(ntl.zones$OBJECTID)
countries_GID_0 <- unique(ntl.zones$GID_0)
# add area km2 #
ntl.zones$sf_area_km2 <- st_area(ntl.zones) / 1000000

setwd(wkdir)

##
## roads
##
countries <-  c( 'Niger', 'Nigeria', 'Chad', 'Cameroon' )

require(lwgeom)
road_lk_chad_shp = paste0(wkdir,'/local/gis_data/018_transportation/roads_panel')

if (!file.exists(paste0(road_lk_chad_shp,'.shp'))){
  roadsafr_sf <- st_read(dsn=roads_dir,layer='rds_remi_afr') %>% 
  sf::st_make_valid()
  ## lwgeom::st_make_valid()

    reg_sf <- st_read(dirname(adm0_lake_chad_shp), basename(adm0_lake_chad_shp))
    GID_0 <- reg_sf$GID_0
    
    roadsafr_prj_sf <- st_transform(roadsafr_sf,crs(reg_sf))
    roads <- subset(roadsafr_prj_sf, COUNTRY %in% countries)
    
  write_sf(roads,dirname(road_lk_chad_shp),basename(road_lk_chad_shp),driver='ESRI Shapefile')
} else {
  roads <- read_sf(dirname(road_lk_chad_shp),layer=basename(road_lk_chad_shp))
}   
        
# slope in percent
# 0.00833
Slope_fn <- file.path(wkdir,"local/gis_data/006_elevation/meanslope/meanslope_lkch.tif")
if(!file.exists(Slope_fn)){
  Slope_globe <- raster(file.path(wkdir,"local/gis_data/006_elevation/meanslope/meanslope"))
  Slope_crop <- crop(Slope_globe, alignExtent(fishnet_grid,Slope_globe,snap='out'))
  extent(Slope_crop)
  extent(fishnet_grid)
  # convert 100
  Slope_radians = Slope_crop / 100
  Slope <- Slope_radians*pi/180
  writeRaster(Slope,Slope_fn)
  rm(Slope_globe)  
} else {
  Slope <- raster(Slope_fn)
}


##
## return ((6 * exp(-3.5 * abs(x + 0.05))) * 0.6)
## waldo in km/h
## 

offpath <- function(x) {
  return ((6 * exp(-3.5 * abs(x*pi/180 + 0.05))) * 0.6)
}

## radians to degrees

onpath <- function(x) {
  return (6 * exp(-3.5 * abs(tan(x) + 0.05) ) )
}

grid_res_km = 11.1

# input: slope in degrees 
# output: km / h 
slope_kph <- raster::calc(Slope, onpath)

# cost in hours per cell
# time = distance / speed
# raster::area no weights
slope_hr.prj <- projectRaster(slope_kph,fishnet_grid, res = 0.1, method='bilinear')
slope_hr.prj <- mask(slope_hr.prj,fishnet_grid)
slope_hr.cell = sqrt(raster::area(slope_hr.prj)) / slope_hr.prj

# by period
roads.t = subset(roads,select=c(COUNTRY,R2008))
# transform road classes into travel cost to cross a pixel
# grid resolution in km

# max speed for each cell 
roads.t$speed = unlist(lapply(roads.t$R2008,FUN=roadcat2speed))
summary(roads.t$speed)
require(stars)
roads.grd.fn = paste0(wkdir,'/local/gis_data/018_transportation/road_speed_all.tif')
if(!file.exists(roads.grd.fn)){
  roads.grd.t = raster::rasterize(as(roads.t,Class="Spatial"),fishnet_grid,field = 'speed', fun='max', background=NA)
  writeRaster(roads.grd.t, roads.grd.fn, driver= "GTiff",overwrite=TRUE)
} else {
  roads.grd.t <- raster(roads.grd.fn)
}

# cost in hours per cell
# time = distance / speed
# raster::area no weights
roads.hr.cell = sqrt(raster::area(roads.grd.t)) / roads.grd.t
plot(roads.hr.cell)

# adjustments for raster alignment and resolution #
extent(slope_hr.cell)
extent(roads.hr.cell)
costgrid.stack <- stack(slope_hr.cell,roads.hr.cell)
# least cost of road path and natural path
cost.grd.t <- stackApply(costgrid.stack, c(1,1), fun='min',na.rm=TRUE) 
# mask country only including full border cell #
cost.grd.t <- mask(cost.grd.t,fishnet_grid)
plot(cost.grd.t)
#cost.grd.t
timecost <- paste0(wkdir,'/local/gis_data/018_transportation/timecostgrid_all.tif')
if (!file.exists(timecost)) {
  writeRaster(cost.grd.t,timecost,driver='GTiff',overwrite=TRUE)
}

## update cost-distance grid with background levels from waldo tobler function
## the values in the input raster are costs. 
## take the arithmetic mean of the cost first and then take the reciprocal of that to get the conductance. 
## ie. The traveller experiences half of the cost of the origin cell and half of the cost of the destination cell, (cost1 + cost2)/2.
## function(x) 1/mean(x)  NOT function(x) mean(1/x)
cost.t = gdistance::transition(cost.grd.t, transitionFunction=function(x) 1/mean(x), directions=8)

##
## END COST
##

##
year <- 2008

  # bhbr
  bhbr_grid <- rasterize(bhbr,fishnet_grid,field="OBJECTID",fun=mean)
  bhbr.df <- as.data.frame(rasterToPoints(bhbr_grid))  
  names(bhbr.df)[3]<- 'bhbr'
  tt_cost.hour <- gdistance::accCost(cost.t, as.matrix(bhbr.df[,1:2]))
  print(paste0("Travel time (open border): ",names(tt_cost.hour)))
  tt_cost.hour <- mask(tt_cost.hour,fishnet_grid)
  writeRaster(tt_cost.hour,paste0(wkdir,'/proc_data/tt2bhbr.tif'),driver='GTiff',overwrite=TRUE)
  ## [tt02bhbr] : travel time (open border) 2 BH between rivers ##
  ntl.zones$tto2bhbr <- exact_extract(tt_cost.hour,ntl.zones,fun='mean')
  summary(ntl.zones$tto2bhbr)
  names(ntl.zones)
  ## add identifiers to grid ##
  tt.df <- subset(as.data.frame(ntl.zones),select=c(OBJECTID,tto2bhbr))
  dim(tt.df)
  save.dta13(tt.df,paste0(wkdir,'/proc_data/tt2bhbr.dta'))


  ##
  ## network based controls
  ## travel time to Ndjamena
  ## travel time to primary city in country (using new def for Cameroon)
  ## travel time to capital city in country (using new def for Cameroon) 
  ##
  
  ## africapolis, poly
  afripolis <- st_read(dsn=dirname(afripolis_shp),layer=basename(afripolis_shp))
  afripolis <- subset(afripolis,ISO %in% c('CMR','NER','NGA','TCD'))
  afripolis_pt = st_centroid(subset(afripolis,select=c("pop2015","pop1950","pop1960","pop1970","pop1980","pop1990","pop2000","pop2010")))
  names(afripolis_pt)
  
  # select capitals from afripolis data #
  cities = c('Abuja','Douala','Lagos','Niamey','Ndjamena/Kouss?ri [TCD]','Yaound?')
  city_control = subset(afripolis,Name %in% cities)
  control_pt = st_centroid(city_control)

## sp join with objectid ##
tt <- stack()
for (i in c(1:dim(control_pt)[1])) {
  print(i)
  loc_i <- control_pt[i,]
  # africapolis city ID
  loc_i_grid <- rasterize(loc_i,fishnet_grid,field="ID",fun=mean)
  loc_i.df <- as.data.frame(rasterToPoints(loc_i_grid))  
  loc_i_name <- as.character(control_pt$Name[i])
  if (loc_i_name=="Ndjamena/Kouss?ri [TCD]") { 
    loc_i_name <- "Ndjamena"
  }
  if (loc_i_name=="Yaound?") { 
    loc_i_name <- "Yaounde"
  }
  
  names(loc_i.df)[3]<- loc_i_name
  tt_cost.hour <- gdistance::accCost(cost.t, as.matrix(loc_i.df[,1:2]))
  print(paste0("Travel time (open border): ",names(tt_cost.hour)))
  #[tt_cost.hour=="Inf"]<- NA
  tt_cost.hour <- mask(tt_cost.hour,fishnet_grid)
  writeRaster(tt_cost.hour,paste0(wkdir,'/proc_data/tt2',loc_i_name,'.tif'),driver='GTiff',overwrite=TRUE)
  names(tt_cost.hour)  <- loc_i_name
  ntl.zones[paste0('tt2',loc_i_name)] <- exact_extract(tt_cost.hour,ntl.zones,fun='mean')
  
  tt <- stack(tt,tt_cost.hour)
}

##
## NEED TO INCLUDE Raster cells on boundary 
##

names(ntl.zones)
tt.df <- as.data.frame(ntl.zones)
tt.df <- tt.df[,c(2,4,31:36)]
names(tt.df)
write.dta(tt.df,paste0(wkdir,'/proc_data/tt2control_cities.dta'))


##
## PLOT LAKE CHAD
##

##
year <- 2008

# World Bank Lake Chad definition
wblakechad_grid <- rasterize(bhbr,fishnet_grid,field="OBJECTID",fun=mean)
##wblkc.df <- as.data.frame(rasterToPoints(wblkc_grid))  
wblkc.df <- as.data.frame(rasterToPoints(wblakechad_grid))  # change
names(wblkc.df)[3]<- 'wblkc'
tt_cost.hour <- gdistance::accCost(cost.t, as.matrix(bhbr.df[,1:2]))
print(paste0("Travel time (open border): ",names(tt_cost.hour)))
#[tt_cost.hour=="Inf"]<- NA
tt_cost.hour <- mask(tt_cost.hour,fishnet_grid)
writeRaster(tt_cost.hour,paste0(wkdir,'/proc_data/tt2wblkc.tif'),driver='GTiff',overwrite=TRUE)
## [tt02bhbr] : travel time (open border) 2 BH between rivers ##
ntl.zones$tto2bhbr <- exact_extract(tt_cost.hour,ntl.zones,fun='mean')
summary(ntl.zones$tto2bhbr)
names(ntl.zones)
## add identifiers to grid ##
tt.df <- subset(as.data.frame(ntl.zones),select=c(OBJECTID,tto2bhbr))
dim(tt.df)
save.dta13(tt.df,paste0(wkdir,'/proc_data/tt2bhbr.dta'))